---
title: "Part2 AL Reliability"
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      warning = FALSE)
```

```{r load libraries, include= FALSE}
library(tidyverse) 
library(janitor) # clean df
library(readxl) # read in excel files
library(irr) # interrater reliability
library(irrCAC) # gwet's AC
library(htmlTable) # table outputs
```

```{r load data, include = FALSE}
## dataframes are set up so each row is possible AL that was rated/denoted by at least on rater. There is a column that donates the specific tooth/the tooth code, what area the AL was marked in (prenatal enamel, postnatal enamel, or dentine), and a column for each individual or team of raters that contains their rating (RL, MAL, or HAL). 

## individual rater data
load("individ_p2.Rdata")

## team rating data
load("team_p2.Rdata")

## weighted lights function created in part 1 analyses
lights_weighted <- readRDS("path_to/lights_weighted.RDS")
```

```{r color scheme, include=FALSE}
# colors for plots 
mycolors <- c("#284654", "#6881A4", "#778980", "#ECB462", "#B49180","#DECFB8", "#e09d7f", "#C46565", "#773863")

# color blind friendly colors 
my_colors2 <- c("black", "#2271B2", "#3D67E9", "#F748A5", "#359B73", "#D55E00", "#E69F00", "#F0E442")
```


```{r data cleaning and set-up, include=FALSE}
# clean var names (data starts in 2nd row of df)
individ_p2 <- individ_p2 %>% clean_names(case = "snake")
team_p2 <- team_p2 %>% clean_names(case = "snake")

## group team members for my reference 
t1 <- c("rater7", "rater2")
t2 <- c("rater8", "rater6")
t3 <- c("rater5","rater3")

## group all raters for overall calculations
rating_cols <- c("rater7", "rater2", "rater8", "rater6", "rater5","rater3")
rating_teams <- c("rater7_rater2", "rater8_rater6", "rater5_rater3")

## create df of "overall" reliability, which collapses MAL and HAL into AL and RL to RL
# copy over data to preserve original ratings 
overall_ind <- individ_p2
overall_team <- team_p2

# collapse MAL/HAL
overall_ind[rating_cols][overall_ind[rating_cols] == "HAL"] <- "AL"
overall_ind[rating_cols][overall_ind[rating_cols] == "MAL"]<- "AL"

overall_team[rating_teams][overall_team[rating_teams] == "HAL"] <- "AL"
overall_team[rating_teams][overall_team[rating_teams] == "MAL"]<- "AL"

## create df for MAL/HAL/RL weighted rating that does NOT include dentine (due to all being HAL)
weighted_ind <- individ_p2 %>% filter(area != "Dentine") 
weighted_team <- team_p2 %>% filter(area != "Dentine")

## make intensity ratings ordered factors 
weighted_ind[rating_cols] <- lapply(weighted_ind[rating_cols], factor, levels = c("RL", "MAL", "HAL"), ordered = TRUE)

weighted_team[rating_teams] <- lapply(weighted_team[rating_teams], factor, levels = c("RL", "MAL", "HAL"), ordered = TRUE)

## prenatal vs postnatal enamel ratings 
# prenatal enamel
prenatal_overall_ind <- overall_ind %>% filter(area == "Prenatal")
prenatal_weight_ind <- weighted_ind %>% filter(area == "Prenatal")

# postnatal enamel
postnatal_overall_team <- overall_team %>% filter(area == "Postnatal")
postnatal_weight_team <- weighted_team %>% filter(area == "Postnatal")
```

```{r overall indiv, include=FALSE}
# check prevalence 
apply(overall_ind[rating_cols],2,table) 

# overall individuals
overall_ind_light <- kappam.light(overall_ind[rating_cols])$value 
overall_ind_gwet <- gwet.ac1.raw(overall_ind[rating_cols])$est$coeff.val 
overall_ind_agreee <- agree(overall_ind[rating_cols])$value 
```

```{r overall members, include=FALSE}
## Team 1 
overall_t1_light <- kappam.light(overall_ind[t1])$value 
overall_t1_gwet <- gwet.ac1.raw(overall_ind[t1])$est$coeff.val
overall_t1_agree <- agree(overall_ind[t1])$value 
weighted_ind %>% select(all_of(t1)) %>% summarise(rl_hal = length(which((rater7 == "MAL" & rater2 == "RL") | (rater7 == "RL" & rater2 == "MAL"))))

## Team 2 
overall_t2_light <- kappam.light(overall_ind[t2])$value 
overall_t2_gwet <- gwet.ac1.raw(overall_ind[t2])$est$coeff.val 
overall_t2_agree <- agree(overall_ind[t2])$value  
weighted_ind %>% select(all_of(t2)) %>% summarise(rl_hal = length(which((rater8 == "MAL" & rater6 == "RL") | (rater8 == "RL" & rater6 == "MAL"))))

## Team 3 
overall_t3_light <- kappam.light(overall_ind[t3])$value 
overall_t3_gwet <- gwet.ac1.raw(overall_ind[t3])$est$coeff.val 
overall_t3_agree <- agree(overall_ind[t3])$value  
weighted_ind %>% select(all_of(t3)) %>% summarise(rl_hal = length(which((rater3 == "MAL" & rater5 == "RL") | (rater3 == "RL" & rater5 == "MAL"))))
```

```{r intensity individ, include=FALSE}
# check prevalence 
apply(weighted_ind[rating_cols],2,table) 

# weighted lights and gwets to penalize intensity ratings 
weight_ind_light <- lights_weighted(raters = rating_cols, 
                                    weighting = c(0,0.5,1), 
                                    ratings = weighted_ind) 

weight_ind_gwet <- gwet.ac1.raw(weighted_ind[rating_cols], weights = "linear")$est$coeff.val 

unweight_ind_gwet <- gwet.ac1.raw(weighted_ind[rating_cols], weights = "unweighted")$est$coeff.val 

weight_ind_agree <- agree(weighted_ind[rating_cols])$value 
```

```{r intensity members, include=FALSE}
## Team 1 ##
weight_t1_light <- lights_weighted(raters = t1, 
                                    weighting = c(0,0.5,1), 
                                    ratings = weighted_ind) 

weight_t1_gwet <- gwet.ac1.raw(weighted_ind[t1], weights = "linear")$est$coeff.val 

unweight_t1_gwet <- gwet.ac1.raw(weighted_ind[t1], weights = "unweighted")$est$coeff.val 

weight_t1_agree <- agree(weighted_ind[t1])$value 


## Team 2 ##
weight_t2_light <- lights_weighted(raters = t2, 
                                    weighting = c(0,0.5,1), 
                                    ratings = weighted_ind) 

weight_t2_gwet <- gwet.ac1.raw(weighted_ind[t2], weights = "linear")$est$coeff.val 

unweight_t2_gwet <- gwet.ac1.raw(weighted_ind[t2], weights = "unweighted")$est$coeff.val 

weight_t2_agree <- agree(weighted_ind[t2])$value  

apply(weighted_ind[t2],2,table)

## Team 3 ##
weight_t3_light <- lights_weighted(raters = t3, 
                                    weighting = c(0,0.5,1), 
                                    ratings = weighted_ind) 

weight_t3_gwet <- gwet.ac1.raw(weighted_ind[t3], weights = "linear")$est$coeff.val 

weight_t3_gwet <- gwet.ac1.raw(weighted_ind[t3], weights = "unweighted")$est$coeff.val 

weight_t3_agree <- agree(weighted_ind[t3])$value  

# unweighted 
unweight_t3_gwet <- gwet.ac1.raw(weighted_ind[t3])$est$coeff.val 
```

```{r overall teams, include=FALSE}
# check prevalence 
apply(overall_team[rating_teams],2,table)

# overall across teams
overall_team_light <- kappam.light(overall_team[rating_teams])$value 

overall_team_gwet <- gwet.ac1.raw(overall_team[rating_teams])$est$coeff.val 

overall_team_agree <- agree(overall_team[rating_teams])$value  
```

```{r intensity teams, include=FALSE}
# check prevalence
apply(weighted_team[rating_teams],2,table) 

# weighted lights and gwets to penalize intensity ratings 
weight_team_light <- lights_weighted(raters = rating_teams, 
                                    weighting = c(0,0.5,1), 
                                    ratings = weighted_team) 
weight_team_gwet <- gwet.ac1.raw(weighted_team[rating_teams], weights = "linear")$est$coeff.val 

weight_team_agree <- agree(weighted_team[rating_teams])$value 

unweight_team_gwet <- gwet.ac1.raw(weighted_team[rating_teams], weights = "unweighted")$est$coeff.val
```

```{r overview table, include = FALSE}
## create df of scores 
results_df <- data.frame(
  group = 
    append(rating_teams, c("across teams", "across individuals")),
  intensity_weighted_gwet = 
    c(weight_t1_gwet, weight_t2_gwet, weight_t3_gwet, weight_team_gwet, weight_ind_gwet),
  intensity_weighted_light = 
    c(weight_t1_light, weight_t2_light, weight_t3_light, weight_team_light, weight_ind_light),
  intensity_agree_percent = 
    c(weight_t1_agree, weight_t2_agree, weight_t3_agree, weight_team_agree, weight_ind_agree),
  overall_gwet = 
    c(overall_t1_gwet, overall_t2_gwet, overall_t3_gwet, overall_team_gwet, overall_ind_gwet),
  overall_light = 
    c(overall_t1_light, overall_t2_light, overall_t3_light, overall_team_light, overall_ind_light),
  overall_agree_percent = 
    c(overall_t1_agree, overall_t2_agree, overall_t3_agree, overall_team_agree, overall_ind_agreee))

## round all values 
results_df[,2:7] <- apply(results_df[,2:7],2,round,2)

## second version without lights kapaa
results_df2 <- results_df[,c(1,2,4,5,7)]

# flip to make each group the colnames (take transpose)
results_df2_flip <- as.data.frame(t(results_df2))
colnames(results_df2_flip) <- results_df2_flip[1,]
results_df2_flip <- results_df2_flip %>% 
  mutate(analysis = rownames(results_df2_flip)) %>%
  relocate(analysis, .before = rater7_rater2)
results_df2_flip <- results_df2_flip[2:5,] 
```

### Overview of Results table
```{r PRINT overview table, include=TRUE, echo=FALSE}
results_df2_flip %>%
  addHtmlTableStyle(col.rgroup = c("none", "#F7F7F7"),
                    col.columns = c("none", "#F1F0FA"),
                    align = "lllll") %>%
  htmlTable(rnames = FALSE,
            caption = "Reliability and Percent Agreement for Identifying Presence of ALs and Intensity of ALs")
```

